RAL-TR-1998-067 
IFUM-633-FT 
[hep-ph/9809451 



Matrix Element Corrections to Parton Shower 
Simulations of Heavy Quark Decay 



G. Corcella 

Dipartimento di Fisica, Universita di Milano and INFN, Sezione di Milano 
Via Celoria 16, 20133 Milano, Italy 

and 

M.H. Seymour 

Rutherford Appleton Laboratory, Chilton, 
Didcot, Oxfordshire. 0X11 OQX. U.K. 



Abstract 

Parton showers are accurate for soft and/or collinear emission, but for a 
good description of the whole of phase space they need to be supplemented 
by matrix element corrections. In this paper, we discuss matrix element 
corrections to the decay t — > Wb and apply our results to the HERWIG 
Monte Carlo event generator. The phenomenological results show marked 
improvement relative to previous versions and agree well with the exact first- 
order matrix-element calculation. 
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1 Introduction 



Heavy quark physics is at present one of the main subjects of theoretical and experi- 
mental particle physics]]] and one of the objectives of future experiments at the LHC 
and NLC. Although total cross sections and other inclusive quantities can be reliably 
calculated using fixed-order perturbative QCD, for analyses that require a detailed de- 
scription of final state properties one needs to sum large logarithmic contributions to all 
orders and include non-perturbative hadronization effects. Monte Carlo event generators 
use parton shower models to perform this resummation and phenomenologically-inspired 
models of the hadronization process. 

The parton shower approach relies on the universality of the QCD matrix elements 
in the dominant regions of phase space: soft and/or collinear emission 0. Multi-parton 
final states are built up by starting from a few well-separated hard partons and se- 
quentially adding one extra parton at a time, distributed according to the soft and 
collinear factorization formulae. In the collinear limit it is straightforward to formulate 
this as a probabilistic Markov chain, as the factorization formulae are directly written 
in terms of probability distributions. In the soft limit, the factorization theorem applies 
at the amplitude level, and many interfering amplitudes contribute equally, so it would 
appear that no probabilistic approach could be formulated. However, the remarkable 
result is that the interference is entirely destructive outside angular-ordered regions 
of phase space. Therefore parton shower algorithms can be extended to correctly treat 
both regions, simply by using opening angle as the ordering variable within the collinear 
approach. This is the basis of coherence- improved parton shower algorithms such as 
HERWIG Q . 

Coherence is also important in setting the initial conditions of the parton shower: 
each hard parton may only radiate into an angular region extending as far as its 'colour 
partner' ||. The assignment of colour partners is controlled by the hard matrix element, 
and thus the properties of the parton shower and the resulting jets are influenced by the 
kinematics of the hard scattering, giving a perturbative origin to coherent phenomena 
like the string effect M. 

It is straightforward to incorporate quark masses into a coherence-improved parton 
shower algorithm J7|. In heavy quark production, the mass terms result in an angle- 
dependent suppression that can be approximated by an angular cutoff: the massive 
quark radiates like a massless one at large angles, and not at all at small angles, 
O^rrig/Eq. Heavy quark decay acts like a new hard process in which gluons are emitted 
coherently by the heavy quark and by its lighter quark decay product. In principle, 
this separation into production and decay phases is only valid for gluon energies above 
the decay width of the quark, which for top quarks is of order 1 GeV, and below that 
the coherence of the whole production and decay chain should be taken into account ||. 
In practice however, 1 GeV is close enough to the cutoff that terminates the parton 
shower that these effects may be neglected in the parton shower and only need to be 
incorporated in the non-perturbative hadronization phase. 

Although the soft and collinear regions dominate the emission probabilities, many 
experimental observables are most sensitive to hard non-collinear emission. A poignant 
example is the kinematic reconstruction of the mass of decaying top quarks, in which the 
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treatment of additional jets can make a sizeable difference, leading to an uncertainty in 
the final result that is almost as large as the experimental uncertainty 0. Unfortunately, 
there is no guarantee that a parton shower algorithm will describe such emission well. It 
should, however, be well-described by the next-order (but still tree-level) matrix element. 
Thus, if one wants to have a good description of the bulk properties of events and the 
rare hard wide-angle emissions simultaneously, one must combine the parton shower and 
matrix element approaches to provide a process-specific 'matrix element correction' to 
the parton shower. 

In general there are two shortcomings to be corrected: it is possible that owing to 
the angular-ordering condition, there are regions of phase space that are not covered at 
all by the parton shower (which we call 'dead zones'); and it is possible that even in the 
region it does cover, the parton shower is not a good approximation to the exact matrix 
element. General solutions to both these problems (which we call the 'hard' and 'soft' 
corrections respectively) were described in Ref. and have been incorporated into 
HERWIG for e + e~ annihilation j|TT| and DIS[0. Similar issues were discussed for the 



JETSET/PYTHIA parton shower algorithm in Ref. [ 13] . 

In this paper we discuss the implementation of matrix element corrections to heavy 
quark decay in HERWIG and show first results. In section 2 we recall the relevant 
features of HERWIG 's parton shower algorithm. In sections 3 and 4 we discuss the 
hard and soft matrix element correction respectively. In section 5 we present some 
phenomenological distributions for top quark decay and compare them with the exact 
QCD results and those obtained from the previous version of HERWIG. Finally in 
section 6 we make some concluding remarks. 



2 The parton shower algorithm 

For the decayQ t — ► bW, the elementary probability to emit an additional gluon is given 
by 

To 

where T is the width for the Born process. In the soft and collinear limits, this can be 
approximated by the universal elementary probability for the emission of an additional 
parton in any hard process: 

iP = d 4 s» P(z)iz • ( 2 ) 

T 2tt A s (q 2 ,q 2 ) 

The HERWIG parton shower is ordered according to the variable q 2 = E 2 £, where E is 
the energy of the emitter, £ = Jr^r- and q\$ being the four-momenta of the produced 
partons; z is the energy fraction of the emitted parton in the showering frame. For 



1 In fact we use exactly the same method for all heavy (bottom and above) quarks, but for clarity 
we always use the language of top quark decay. Note also that although we describe the decay as 
t — > bW(g), we actually use the full three-body matrix elements for t — > bff' so that the //' mass 
distribution and angular correlations are correct. 
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massless partons we have £ = 1 — cos 6, where 9 is the opening angle of emission, so that 
ordering in q 2 corresponds to ordering in angles. 

The crucial quantity A s (q 2 , q 2 ) is the Sudakov form factor, representing the probabil- 
ity that no resolvable radiation is emitted from a parton whose upper limit on emission 
is q 2 . The resolution requirement, which in the case of HERWIG is a cutoff on transverse 
momentum, is q 2 . Thus the ratio of form factors appearing in Eq. (fj) represents the 
probability that the emission considered is the first i.e. the highest in q 2 . In diagram- 
matic language, this actually sums up the contributions of virtual and unresolvable real 
emission to all orders. 

Using HERWIG's definitions for q 2 and z, the results of the parton shower are not 
Lorentz covariant, although they become so in the exactly soft or collinear limits in 
which the algorithm is formally valid. Nevertheless, the extrapolations away from those 
limits depend on the frame in which the shower is developed. HERWIG uses a different 
frame for each jet i, in which £j max = 1 and hence, since q 2 = E 2 £, E { = gj max . After 
parton showering, the momenta of the jets cannot be identical to the momenta of the 
partons that initiated them, since the latter are on mass-shell, so some momentum must 
be transferred between jets to conserved energy-momentum overall. The precise details 
become irrelevant in the soft and collinear limits in which the algorithm is valid, but are 
again important for the extrapolation away from that limit that we are concerned with 
here. The procedure is again frame-dependent and is performed by Lorentz boosting 
each jet along its direction in the rest-frame of the hard process that produced it. 

Colour coherence in the hard process dictates the value of q^^ which for colour- 
connected partons i and j are related by gVaxQjmax = Pi • Pj- In principle, this is the 
only requirement, and one has a free choice for the individual values gj max and gj max , 
corresponding to different choices of frame, but in most cases, symmetric limits are 
chosen: qf max = gj max = Pi • Pj- However, for the special case of top quark decay, we 
choose E t = qWax = tth and therefore Ef, = <2W ax = E ^ k , corresponding to the top 
rest-frame, in which the top quark itself does not radiate. 

Because of this choice, and the fact that HERWIG only radiates into the angular- 
ordered region £ < 1, there is no radiation in the W hemisphere, while in the full matrix 
element such radiation is suppressed but not absent. Therefore the dead zone actually 
includes part of the soft singularity, in contrast to the e + e~ and DIS cases in which 
soft emission from one or other jets covers the whole of the angular range. This poses 
some additional problems, which we describe heuristically here, and in more detail in 
the following sections. 

Since HERWIG's dead zone includes part of the soft singularity, the total amount of 
emission into it, if calculated naively, is infinite. Therefore the hard correction method 
of Ref. jnj cannot be used. The correct solution to this problem is to modify HERWIG's 
parton shower algorithm so that it populates the backward direction. Unfortunately, 
this is far from straightforward to do. Instead, we look for simpler, more approximate, 
solutions. 

If we apply a cutoff on the gluon energy then, provided the total probability of 
emitting into the dead zone above the cutoff is small, the hard correction can be used. 
Although this is only an approximation we can easily check, by varying the cutoff, how 
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good it is. In section 5 we show that varying the cutoff in a reasonable range gives a 
negligible effect on physical distributions. 

However, the fact that we apply an energy cut in the hard correction has implications 
for the soft correction. The general idea of HERWIG's approach to modelling emission 
from heavy quarks |]] is that in the soft limit, the emission pattern is similar to that 
from a light quark at large angles, but is smoothly suppressed at small angles, according 
to the replacement 

1 , 02 (3 ) 

where 7 is the Lorentz factor of the emitting quark. This is approximated by a step 
function at 6 = I/7. While this does not give a very good description of the angular 
distribution of emitted gluons, it does give a good approximation to the total amount 
of gluon emission, which is important for the total energy loss of the heavy quark for 
example. In other words, we get the amount of gluon emission about right, but tend to 
put it in the wrong place. For top decay, this is at its most extreme, since we work in 
its rest-frame in which the neglected angular region is the whole of the W hemisphere. 

If we were to apply a soft matrix element correction in the usual way, we would 
decrease the amount of emission in the b hemisphere but, because we use a cutoff on 
the gluon energy in the hard correction, for gluons below this cutoff there would be no 
corresponding increase in the W hemisphere. We would therefore get the total amount 
of emission wrong, and hence quantities like the b quark's energy spectrum would be 
poorly described. To avoid this, we use the same cutoff in the soft correction, and only 
correct the distribution of gluons with energy above it. Below the cutoff, the standard 
parton shower description is used alone. 

In summary, we can construct the matrix element corrections as the combination of 
hard and soft corrections in exactly the usual way, but only in the region of gluon energies 
above a given cutoff, E min . Below the cutoff, we use the parton shower uncorrected. We 
use 2 GeV as the default value of the cutoff. 

Finally, as discussed in Ref. ||10|| , one should include a form factor in the region 



covered by the hard correction, to ensure smooth matching at the boundary. Since 
we approach the soft singularity in the top decay case, this could in principle be more 
important than in other cases. However, we find that even with a cutoff as small as 
1 GeV, the form factor is never more than a few percent below unity, so does not 
significantly affect any physical distributions. We therefore neglect it for simplicity, as 
is already done in the e + e~ and DIS cases. 



3 Hard Corrections 

The main step in constructing both the hard and soft matrix element corrections is 
to relate the variables generated by HERWIG, z and £, to the kinematic variables we 
use in the matrix element calculation. We parametrize the phase space of the process 
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t(q) — > W{j>i)b(p2)g(p-s) in terms of the variables^ 



2p 2 -p 3 2pt-q 

Xl = l 2~ = 2 ° ' ( 4 ) 



where we have defined^ a = mf^/m 2 , and 



x 3 = — • ( 5 ) 

The phase space limits for the decay t — > oVTg are as follows: 

-^_ + (l-z 3 ) < X! < 1, (6) 

1 — x 3 

< :r 3 < 1 - a. (7) 

In order to express x\ and 23 in terms of £ and z, we observe that the Lorentz boost of 
the b quark's parton shower takes place in the top rest frame, so that the mass, m, of 
the b—g jet and its transverse momentum, q t , relative to the b-W axis are conserved. We 
also define the energy of the b quark in the showering frame to be y/cmt with, because 
HERWIG uses the top rest frame, c = (1 — a) 2 /4. In terms of the showering variables, 
we obtain: 

m 2 = 2z(l-z)f cm 2 , (8) 

and in terms of the matrix element variables: 

m 2 = m 2 (l-xi), (10) 

2 m 2 (1 ~ glKgl + ^3(2 - - a) - x| - 1) 

% * (xi + a) 2 -4a ' 1 J 

They can be combined to give: 



1 r J( Xl - (1 -c))/c 1 

^ 1 + V (2x 3 + an - 2 + a) , (12) 

; (x! + a) 2 - 4a J 



, = 2(l-x 1 )((x 1 + a) 2 -4a) 

" (1 - xi)(2x 3 + Xi - 2 + a) 2 + 4c(l - x 3 )(x 1 + x 3 - 1) - 4acx 3 ' 

The region not populated by HERWIG can be calculated from the condition £ > 1: 

< x 3 < 1 - a, (14) 

(1 - £3) + ^ 3 < £1 < £l max , (15) 



2 To simplify the formulae in this paper, we neglect the b mass, although we do include it in HERWIG, 
and hence in all figures that we show. 

3 Recalling that we use the three-body matrix elements for t — > bff, a is actually defined to be 
(Pf + Pf') 2 / m ti which varies from event to event, but this difference is unimportant here. 
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Figure 1: Phase space limits on x\ 



— — 



a and x 3 



2P3Q 



in the decay t — > bWg 



(solid) and the edge of the region covered by HERWIG (dashed). Also shown, barely 
visible in the top-left corner, is a 2 GeV cutoff on the gluon energy. 



with Xi max being the solution of the cubic equation 

x\ + (3 + 2a - Ax 3 )xj + (a 2 - 4(1 - x 3 )(2 - c- x 3 ) - 2a(3 + 2x 3 ))x x 

+a(4 - a) (2 - c) + (1 - c)(2(l - x 3 ) - a) 2 

Such a solution can be expressed in the following form: 



x u 



23ft < r + iJq 3 - r 2 



1/3) 1 

--(3 + 2a-4x 3 ) 



0. (16) 



(17) 



where we have defined 



q= l -CZ + 2a-Ax 3 ) 2 - 1 - 



a 2 - 4(1 -x 3 ) (2 — c — x 3 ) — 2a(3 + 2x 3 ) 



and 



(18) 



r = -(3 + 2a-4x 3 ) [a 2 - 4(1- x 3 ) (2 - c - x 3 ) - 2a(3 + 2x 3 ) 

- l - (a(4 - a) (2 - c) + (1 - c) [2(1 - x 3 ) - a} 2 ) - i(3 + 2a - 4X3) 3 . (19) 

In Fig. we plot the full phase space limits and HERWIG 's limits. Note that the collinear 
limit is X\ — 1 and the soft limit is x 3 = and that, as anticipated, HERWIG's missing 
region includes part of the soft limit, which we avoid by using the cutoff x 3 > 2E m i n /m t . 

In the dead zone we use the exact differential width, which, continuing to neglect 
the b mass, is: 



1 a^r 

r dx\ dx 3 



1 



a 



l-o)(l + i-2a 



1 



xUl - x x . 



(l+i-2o)[(l-o)a!3-(l-a!x)(l-a;3)-a!i] 



+(l+ir-)x 3 ,(x 1 +x- i -l) 2 +2xl{l-x 1 ) >, 



(20) 



It is straightforward to use standard techniques to generate this distribution within the 
allowed phase space region. 

A final point involves the scale at which we evaluate a s . As is well known, using the 
transverse momentum of the emitted gluon sums up an important class of next-to-leading 
logarithmic corrections ]14]|. However, we should be careful in using the transverse mo- 



mentum in the backward hemisphere, because in the W direction it is zero, even though 
the gluon is not collinear to any coloured parton. We instead use the Durham-like 
generalization of transverse momentum, 

^ = 2E 2 (l - cos**) = - X ^~ X ^ ml (21) 
y 2 — x\ — £3 — a 

4 Soft Corrections 

According to Ref. [O, the soft correction is implemented by multiplying the emission 
probability for every emission that is the hardest so far by a correction factor that 
is simply the ratio of HERWIG's differential distribution to the matrix element one. 
The only non-trivial part of this is in calculating the Jacobian factor J(xi,x^) of the 
transformation (z, £) — > (x\,x$). HERWIG's distribution is then given by 

d 2 T _ d 2 T ^ 
dxidxz dzdS, 

where d 2 T jdzd^ is given by the elementary emission probability given in Eq. (H). J can 
be simply calculated from the relations given earlier, but we have not found a particularly 
compact form so we do not reproduce it here. 

We recall that we only apply the soft correction for x 3 > 2E min /m t . 



5 Results 



The exact matrix element for e + e — > (bW + ) (bW ) g via top quarks was calculated 
in Refs. |i5| , |l6H and compared to HERWIG versions 5.8 and 5.9 respectively. In fact, 



serious bugs have been found in the treatment of top decays in both versions, so no firm 
conclusion could be reached about the need for matrix element corrections. 

In this section, we repeat the analysis of Ref. [16] in which they consider a centre- 
of-mass energy only slightly above threshold, = 360 GeV, so that essentially all 
gluon emission is associated with the top quark decays. We work at the parton level 
(the final state of the parton shower, without hadronization or b decay) and cluster all 
partons into three jets using the k± algorithm]!?]]. Then the jets are required to be hard, 
Exj > 10 GeV, and well-separatedf], AR > 0.7. We force the W decays to be leptonic 
and do not include their decay products in the jet clustering. 

In Figs. ^| and we show the differential distribution with respect to AR, the 
separation of the closest pair of jets in the event, and log 2/3, where y 3 is the value of the 

4 Note that the text of Ref. says that the angular cut is AR > 0.4, but it is clear from their 
Fig. 7(a) that it is actually AR > 0.7. 
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AR 

Figure 2: Differential distribution of in- 
variant opening angle, AR 2 = Arf + Acf 2 
for three-jet e + e~ — > ti events at y/s — 
360 GeV, according to HERWIG version 
5.9 (dashed), 6.0 (dot-dashed) and the 
version with matrix element corrections, 
6.1 (solid). 



i°g(y 3 ) 

Figure 3: As Fig. |2| but for the distri- 
bution of Durham jet algorithm measure, 
y 3 = min^ff mm(E 2 , E 2 )(l - cos%)]. 



cutoff in the k± algorithm at which the three jets would be merged to two, according to 
the last public version, 5.9, version 6.0, which has the bug that was found in 5.9 fixed, 
and a preliminary version that includes our matrix element corrections, 6.1. The version 
5.9, which was used in Ref. [116] , clearly gives too much gluon radiation. This is easily 
explained by the bug which was corrected in version 6.00, which is seen to have much 
less radiation at large angles and y 3 values. 

The only difference for top decays between versions 6.0 and 6.1 is in the matrix 
element correction discussed here. Therefore Figs. ^| and ^| directly show its effect. 
At large angle and y 3 the rate of gluon emission is greatly enhanced due to the hard 
correction, which fills the missing region of phase space. At smaller angles and y 3 values, 
the soft correction gives a small reduction, showing that in the uncorrected version the 
deficit at large angles was at least partially compensated by a surfeit at small angles. 

We turn now to a comparison with the tree-level matrix element results of Ref. [|ll| . 
In order to remove dependence on the electroweak production process (for example the 
fact that HERWIG uses the cross section for on-shell top quarks while the matrix element 
calculation implicitly includes the top width) we compare the differential distribution 
normalized to the total e + e~ — > ti cross section. Even after doing this, the normalization 
of a tree-level calculation is still rather arbitrary, since the scale of a s is not fixed. We 
choose to set a s equal to its value at a typical transverse momentum of jets passing 
the cuts, a s (~ 30 GeV) ~ 0.145. Results are shown in Figs. | and |^. After applying 
matrix corrections the agreement is very good at large y 3 and AR. As y 3 is reduced, we 
probe lower momentum scales, where a s is getting larger in HERWIG, while it is fixed 
in the matrix element calculation, so HERWIG lies above the calculation. Finally, at 
very small y 3 , Sudakov suppression from multiple emission would be expected, which is 



5 The wrong frame had been used for the b parton shower: heavy quark decay was treated like all 
other processes and showered in the frame in which Eb = ^Jp t ■ Pb ~ 110 GeV instead of the top rest 
frame in which E b = » 70 GeV. 
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AR 

Figure 4: As Fig. g but from HERWIG 
6.1 (solid) and the tree-level calculation 
with a s = 0.145 (dashed). 



-4 -3.5 -3 -2.5 -2 -1.5 -1 

i°g(y 3 ) 

Figure 5: As Fig. | but from HERWIG 
6.1 (solid) and the tree-level calculation 
with a s = 0.145 (dashed). 



included in HERWIG but not in the matrix element calculation, so HERWIG lies below 
the calculation. 

Finally, in Figs. ^ and [7| we show the dependence on the arbitrary cutoff on soft gluon 
energies we introduced in the hard matrix element correction. It is clearly insignificant. 
It is worth noting that the actual fraction of events that have an emission in the dead 
zone generated by the hard matrix element correction varies considerably with the cutoff, 
from 2.4% at 5 GeV to 3.8% at 1 GeV, so the fact that the physical distribution after 
cuts is unaffected is a non-trivial test of the self-consistency of the procedure. The 
default cutoff value, 2 GeV, was used for the previous plots. 




Figure 6: As Fig. | but from HERWIG 
6.1 with gluon energy cutoff set to 1 GeV 
(dashed), 2 GeV (solid) and 5 GeV (dot- 
dashed) . 



Figure 7: As Fig. | but from HERWIG 
6.1 with gluon energy cutoff set to 1 GeV 
(dashed), 2 GeV (solid) and 5 GeV (dot- 
dashed) . 
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6 Conclusions 



We have discussed HERWIG's parton shower treatment of top quark decay. We applied 
the exact matrix-element result for emission in the region of phase space that it does 
not populate and for the hardest emission in the usually-populated region. Since the 
missing region includes part of the soft singularity, we have had to apply an arbitrary 
cutoff on the energy of gluons emitted in this region. The dependence on this cutoff is 
however negligible after applying jet cuts. 

HERWIG's predictions are considerably improved by this correction, and agree well 
with the exact leading-order perturbative results. Where there are differences, they can 
be understood as advantages of the parton shower approach. We therefore feel confident 
that HERWIG now gives a reliable treatment of top quark decays. 

It is clearly interesting to ask how much these corrections affect the top mass re- 
construction from final state properties at the Tevatron. However, before doing so it is 
essential to also implement matrix element corrections to the top production processes. 
This work is in progress. 
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